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This  thesis  implements  an  ‘‘optimar  Wald  sequential  hypothesis  testing  scheme 
for  the  general  diffusion  model, 


dx,  =  f'{Xt)dt  +  g'{xi)dV/\  (1.1) 

dy,  =  /i'(x,)dt  +  7dV,  (1.2) 

where  i  =  0, 1.  Under  hypothesis  i; 


X,  is  an  n‘  dimensional  state  vector; 

Wj  is  an  m'  dimensional  Brownian  motion  vector; 
y,  is  a  p  dimensional  observation  vector; 

V,  is  a  p  dimensional  Brownian  motion  vector; 
identical  under  both  hypothesis; 

f',  g',  h'  are  known  functions  of  the  respective  state  vectors. 

So  f'{Xt)  is  an  n’  vector,  g'{Xt)  is  an  n’  x  m‘  matrix,  and  /i*(xj)  is  a  p  vector. 
Equation  1.1  is  called  the  state  equation  and  equation  1.2  is  the  observation 
equation.  When  necessary,  we  will  drop  the  time  subscript  and  use  subscripts 
to  denote  elements  in  vectors  or  matrices. 

Given  observation  data  from  one  of  the  hypotheses,  we  wish  to  determine 
from  which  hypothesis  the  data  came.  Quoting  results  from  [LS78)  we  will 
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present  an  "optimal”  sequential  detection  scheme  and  then  formulate  a  solution 
with  several  possible  numerical  implementations  of  the  detector.  We  will  give 
some  numerical  results  and  compare  several  of  these  numerical  schemes  through 
computer  simulation. 

To  begin,  let  us  review  some  results  from  chapter  17  of  [LS7S].  Consider 
the  following  pair  of  stochastic  differential  equations 

Ho  dyt  =  dV, 

Hi  dyt  =  h{xt)dt  +  dV,  !/o  =  0 

These  correspond  to  scalar  observation  equations  for  a  single  hypothesis  versus 
noise  case.  Let  be  a  given  probability  space  with  a  nondecreasing 

family  of  cr-algebras  Pt  where  t  >  0  with  P’1  C  P.  Furthermore,  let  W  = 
{Wt,Pt)  be  a  Wiener  process.  Assume  h{xt)  =  {h{xi),Pt)  is  an  unobservable 
process  independent  of  W.  We  make  the  following  definitions  which  describe 
the  sequential  detection  scheme. 

Definition  ;  A  sequential  hypothesis  testing  scheme  is  denoted  b y  A  =  A(r,  S). 
T{y)  signifies  the  decision  time  and  6{y)  6  {0, 1}  the  decision  with  6  =  1  signify¬ 
ing  the  oijservcd  data  came  from  the  model  corresponding  to  hypothesis  1  and 
(>  =  0  from  hypothesis  0. 

Definition  :  We  define  an  error  of  the  first  kind  as  q(A)  =  Py{6{y)  =  0) 
which  in  radar  terminology  is  called  the  probability  of  miss  and  denoted  by  P,\/. 
We  define  an  error  of  the  second  kind  as  P{A)  =  Po(^(y)  =  1)  which  in  radar 
terminology  this  is  called  the  probability  of  false  alarm  and  denoted  by  Pf. 

Definition  :  We  define  Aq,^  as  the  class  of  schemes  A  =  A(t,  6)  with  a(  A)  <  Q 
and  /?{A)  <  p  where  a,p  are  constants  with  a  -f  /?  <  1. 

Definition  ;  Define  hi  =  E\  ih{xt)\pf) 

We  make  the  following  assumptions  required  for  the  proof  of  the  following 
theorem, 
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Al  £,  [|/i(.r()l]  <  oo  f  <  oo  .  i  =  0. 1 

A2  P,  {/o’®  ii'ids  =  oc}  =  1  i  =  0, 1 

A3  Ei  [/o'  /Vff/sj  <  oo  t  <  oo  ,  j  =  0, 1 


Theorem  1  Given  assumptions  Al  to  A3  then  in  the  class  Aqj  there  exists 
a  scheme  A  =  A{f,6)  optimal  in  the  sense  that  for  any  other  scheme  A  = 


A(r,6)  e  Aa,^ 


z  =  0,l  . 


The  scheme  A  =  A(f, 6)can  be  defined  by  the  relationship 


f{y)  —  infimum  {t  >  0  :  A,(t/)  0(.4,  B)} 


where 


and 


Hy)  = 


0  Af(y)(y)  >  B 
1  Af(y)(y)  <  A 


At(y)  =  exp  hsdys  -  ^ 


1-/3  ’  P 

At  is  known  as  the  likelihood  ratio  and  ln(A()  is  the  log  likelihood  ratio.  In  this 
case 


[L 


dt 


<00  i  =  0, 1 


Lemma  1  For  the  scheme  A  =  A(f,  ^), 


Wiy)  <oo)  =  Fi(f(7j)  <  oo)  =  1 
and  a(A)  =  a, /3(A)  =  p. 


The  proofs  for  the  above  theorem  and  lemma  are  given  in  Chapter  17  of  [LS78] 
for  the  scalar  case. 

The  proofs  for  the  equivalent  theorem  and  lemma  for  the  vector  case  are 
given  in  [LaV86|.  For  the  vector  case  assumptions  Al  through  A3  become: 
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Bl  £,  lll/i(X()l|]  <00  t  <  00  ,  i  =  0, 1 
B2  P,{!^\\h,\\^ds  =  oc}  =  l  i  =  0,l 
B3  Ei  [/o  ll/i,|l^ds|  <  oo  t  <(x>  ,  i  =  0, 1 
and  then 

A,(y)  =  exp  jjf  h^dy,  -  |  ^  Ilh,||-d5  j 

where  ^  in  hj  denotes  vector  transpose. 

For  the  general  hypothesis  testing  problem  between  two  hypotheses,  corre¬ 
sponding  to  eqs.  1.1  and  1.2,  the  likelihood  ratio  is  given  by 

.\t{y)  =  cxpj^  (h\  -  h^f  dy,  -  [\\h\\\^  -  ||A°||2)  dsj  . 

With  the  additional  assumption 

Pi  \\iA  -  =  ooj  =  1  i  =  0, 1 

we  guarantee  a  finite  decision  time. 

In  order  to  calculate  A((y)  it  is  necessary  to  determine  h\,  the  conditional 
expectation  of  /i’(x<)  given  the  observation  y,  s  <  t.  Under  appropriate  condi¬ 
tions,  the  conditional  density  of  Xj  conditioned  on  the  observation  dy,,  which  we 
will  denote  as  «,(x,  t),  satisfies  the  linear  parabolic  partial  differential  equation 
known  as  the  Zaliai  equation  (BBH83].  In  the  next  chapter  we  present  the  Zalcai 
equation  and  the  necessary  conditions  for  existence  of  a  unique  solution  which 
in  turn  will  yield  At{y).  We  also  derive  several  approximation  schemes  which 
are  implemented  and  we  make  comparisons  of  the  results. 
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Chapter  2 


Solution  and  Numerical 
Approximation  of  the  Zakai 
Equation 


2.1  Solution  of  the  Zakai  equation 


The  Zakai  equation  is  given  by  (  2.1).  When  7  (the  observation  noise  scaling 
factor)  is  equal  to  one  the  conditional  density  u{x,t)  of  the  underlying  state  x 
conditioned  on  the  observations  dy  is  known  to  satisfy  the  Zaltai  equation.  We 
will  drop  the  hypothesis  superscripts  and  the  time  subscripts,  instead  using  sub¬ 
scripts  to  denote  vector  and  matrix  elements.  So  the  functions  /i(x,),  /q,  y/,  and 
Xt  will  not  have  explicit  time  dependence.  The  following  assumptions  guarantee 
existence  and  uniqueness  of  the  solution  to  the  Zakai  equation  [BBH83]  and  are 
assumed  to  hold  throughout  the  chapter. 

Cl  L'  is  uniformly  elliptic 


C2  /(x),ct(x),/i(x),  ^/.(x),^<7,j(x),^(7,j(x),5^a,j(x),  2^/it(x) ,  and 
for  i,j  =  1, ...,  n  and  k  =  1,  ...,pare  uniformly  bounded  and  Lipschitz 

continuous. 
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Then  the  unormalized  conditional  density  u(x,t)satisfies  (  2.1), 


du{x,t)  =  L’u{x,t)dt  -  ^\\h\\^u{x,t)dt  +  li^dyu{x,t)  (2.1) 
L'uix.t)  =  [g..j(x)u(x,t)] 

=  ^<7(x)/(x) 

where  (Ti^{x)  is  the  i,  j  element  of  the  the  matrix  a(x).  To  determine  the  partial 
differential  equation  that  u(x,t)satisfies  for  7  ^  1  it  is  simply  neccesary  to 
renormalize  the  output  equation,  so  h{x)  becomes  ^/j(x)  and  the  observation 
dy  becomes  f^/y. 

We  wish  to  solve  (  2.1)  for  n(x,t).  We  factor  the  term  L‘u{x,  t)  as  A‘u{x,  t)  + 
C{x)u{x,t)  where  ,4*  contains  the  terms  which  have  derivatives  of  the  density 
u(x,t).  This  is  done  based  on  numerical  implementation  considerations  which 
will  be  expanded  upon  later.  Then  (  2.1)  can  be  written  as 

duix,t)  =  .4*u(x,t)dt  +  |c(x)  -  ^l|h|l^j  u(x,t)df  4- /i^dyu(x,t)  .(2.2) 


Defining 

d0(x,y,t)  =  li^dy  +  \ 
we  will  use  a  Gauge  Transformation  with 


CM  -  i||;.ir' 


dt 


r(x,t)  =  exp"'*^*’^’'^u(x,t) 


so 

u(x,  t )  =  exp*^*’^’*^  r(x,  t)  . 

Now  n{x,t)  muot  satisfy  (  2.2)  so 

dr(x,  t )  exp^**’^''^  =  A’r{x,  t)  exp'^*’^’*^  dt  +  (dd>(x,  y,  t))  r(x,  t)  e.xp^*’''^’'^ 


exp^*’^’*^  =  dr{x,  t)  +  r(x,  t)d  exp'^*’^’*^  =  w4*r(x,  t)  exp^*’^’*^  dt 


+  {d<f>{x,y ,t))  r{x,  t)  exp'^^*'^''’ 


(lr{xj)  =  exp~‘^*-y'‘’  [  .4’r(x,0exp^*'y’'’  dt  +  {d(p{x,y,t))  r(x.  f) 

-r(x,f)(dexp^*-y''))] 

dr(xj)  =  exp"'^’'  ^''^  [  exp^^*’^’'^  dt  +  {d4>[x,y,t))  r(x, f)  exp'’*’'-^  '’ 

-  {d[(i>{x,  y,  0)  r{x,  t)  exp^^’'-^'"  ] 

dr(x,t)  =  exp"*^*’^’*^  ^>l’r{x,f)  exp*^*’^’**]  dt 
or 

<ir(x,  t)  =  exp""^*’^’*^  -4’  exp^*'^’'^  r(x,  t)dt  . 

This  is  a  classic  parabolic  partial  differential  equation.  Using  results  from  semi¬ 
group  theory  and  differential  equations  the  solution  is  given  by 

r(x,  t)  =  exp“^*’^’'^  exp  '**  exp^*'^’*^  u(x,  0) 

so 

u(x,0  =  exp'^*’^’*^  r(x,t)  =  exp-^**exp'’^*’^’‘’ u(x,0) 

=  exp”^** u(x, 0)  .  (2.3) 

2.2  Calculation  of  Likelihood  Ratio  A 

From  the  introduction  we  have  for  the  single  hypothesis  versus  noise  case 

A,(y)  =  exp  I  hfdy,  -  ^ 

and  for  the  general  hypothesis  case 

A,(y)  =  exp  {hi  -  h^Y  dy,  -  (l|/iill'^  -  \\h°\\^)  dsj  . 

To  calculate  A(  from  the  above  formulas  it  is  necessary  to  find  /i|,  the  conditional 
expectation  of  /i|.  We  will  solve  the  Zakai  equation  for  Ui{x,t)  the  unnormalized 
conditional  density  assuming  hypothesis  i,  then  taking  the  expectation  of  h\  with 
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respect  to  the  normalized  conditional  density  we  can  calculate  A(.  However,  we 
will  show  that 

A((y)=  /  ui{x,t)dx 
Jnv' 

in  the  single  hypothesis  case  and 


Sst,Uo(x,t)dx 

in  the  general  hypothesis  case.  By  using  the  above  formulas  to  calculate  A<  we 
avoid  the  additional  approximation  in  calculating  h\. 

Lemma  2 


A,(y)  =  cxp|^  fi^dy, 

(for  the  single  hypothesis  versus  noise  case)  satisfies  the  stochastic  differential 
equation 

dA,{y)  =  A,{y)h^dy  . 

proof:  Changing  notation  let  At(y)  be  denoted  by  A(y,  t)-  ^se  '  to  denote 
differentiation  and  the  corresponding  subscripts  will  denote  the  variable  with 
respect  to  which  the  differentiation  is  performed.  Let 

:(t}^  fi,:dy.-'-l)\\h£ds 

then 

d2it)  =  h^dy -^\\h\\'^dt  . 

Using  the  fact 


dy  =  h{x)dt  +  dY 
=  h{x)dt  +  dV 

where  c/V  is  a  Brownian  motion  vector  independant  of  dV  but  with  the  same 
distribution  and  is  obtained  by  taking  expectations,  we  get 

=  5I W*  +  . 
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So  z{t)satisfies  a  stochastic  differentail  equation.  Applying  Ito  s  formula  to 

A(2(f),0  =cxp{r(0} 

we  get 

dMz{t)j)  =  +  + 

+/i^A;(2(t),f)rfw 

=  +  51HPA(2(0,<)1  dt  +  /JAiz{t),t)dW 

=  ll/ill2A(-(t),i)dt  +  h'^A{z{t),t)cm 
=  A(2(t),t)  [ll^lpdf  + 

=  A{z{t),t)^hjhdt +  h^c[W^ 

=  Ai z{t),  t)hT  ^hdt  +  dW^ 

Substituting  dy  back  in  for  hdt  +  <iW  we  get 

dA{y,t)  =  A{y,t)h^dy  . 


Consider  the  inner  product  defined  by 


y{x)Ji{x)>=[  gix)h{x)dx 
Jnv' 


For  the  Zaliai  equation  we  have 
L-u(x,t)  =  Yi 

,  OXiOXj  ,_1  OXf 

with  adjoint  operator  L  given  by 

i.(x,<)  =  t  wx,<)i  - ■ 

So  <  L’g{x)Ji{x)  >=<  g{x),Lli{x)  >. 

From  (  2.1)  we  have 

d<  u(x,f),l  >=<  L’u(x,t),l  >  dt+  <  u{x,t)h^,l  >  dy  . 
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Then 


d  <  u(x,  t),  1  > 


<  u{x,t),Ll  >  dt+  <  u{x,t)h^,  1  >  dy 

<  u(x,t),hT  > 


<  u{x,t),l  > 

<  u(x,f),  1  >  h^dy 


<  u{x,t),  1  >  dy 


with  <  u(x,0),l  >=  1.  So  A,  and  <  u(x,t),l  >  satisfy  the  same  stochastic 
differential  equation  and  are  equal  P-a.s.  So 

At  =  <u(x,t),  1> 

=  /  u{x,t)dx  . 

For  the  gcnaral  hypothesis  case  we  get 

^  <ni(x,<),l> 

<  no(x,t),  1  > 

_  /jB,  ui{x,t)dx 
Jjfp,uo(x,t)dx 


2.3  Numerical  approximation  of  the  Zakai  equa¬ 
tion 


Discretizing  (  2.3)  in  time  tvs 


L^J 


U{x,t)  =  exp"'*^* =  u(x, 0) 


(2.4) 


/=o 


LaVignashows  u{x,t).  Similarly,  he  shows  the  approximation 


l^J 


U{x,t)  =  u(x,0) 


(2.5) 


1=0 


also  converges  to  u(x,  t)  as  Af  — *  0. 

We  will  discretize  (  2.4)  and  (  2.5)  in  space  (ie.  with  respect  to  the  underlying 

state  x).  We  will  replace  the  differential  operator  A'  by  matrix  operations. 

We  will  use  boldface  capital  letters  to  denote  matricies.  We  discretize  the 
density  u(x,t)  as  a  matrix 
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u(Xo,f)  u{Xo  +  AX2,t)  ■■■  u{Xo  +  kj\X2,t) 

a(xo  +  Axi,i)  m(xo  +  Axi  +  Ax2,t)  •••  u(xo  +  Axi  +  A.AX2,  t) 

u(xo  +  jA.'ii.t)  ii(xo +  iAxi  +  Ax2,t)  •••  u(xo  +  jAxi  +  AAx2,t) 
which  we  will  denote  b\'  U{x,  t)  where  Xq  is  the  left  endpoint  of  the  discretization 

space. 

Consider  the  first  and  second  order  approximations 

d  ,  u{i  + Ax,,t) -u(x,f) 

— ti(x,t)  = - - - 

ax,  Axi 

d  ,  ^  uix,t)  -  ii{x  +  AXi.t) 

= - - - 

d'^  !  ^  _  u{x  +  AXi,t)  -  2u{x,t)  +  uix  —  AXi,t) 

51?“'*’  ’  A4 

where  Ax,  is  the  discretization  increment  of  the  component  of  the  state  vector 

and  is  scalar  and  AXj  is  an  n  dimension  vector  with  Axi  in  the  position  and 
zero  elsewhere.  So  ^u(x,t)and  ^u(x,t)are  the  second  and  first  derivatives  of 
the  density  with  respect  to  the  t"'  component  of  the  state.  Define 

-2Af  .  .  At  r  /  X  o  ^  /  N 

-  5;-  -/■(*) + 

UDi(x)=  ^  Biax  0,  -/^(x)  +  2  A,,.  j(x) 

LDi(x)  =  ~  ^  min  0,  -/.(x)  +  2^a,,,{x) 

where  the  functions  max[-,-]  and  min[-,-|  take  the  maximum  or  minimum  of 

their  respective  arguments.  We  define  the  tridiagonal  matrices  Aj  and  Aj 
respectively  by 

£>,(Xo)  -Dm(Xo) 

ZDi(Xo+AX,)  D,(Xo  +  aXi)  i;£)i(Xo  +  AXi) 

L£)i(Xo  +  (>-l)AX,)  D,(Xo +  (j  -  l)AXi)  yD,(Xo  +  (j  -  l)AX, ) 
-D,(Xo+jAX,)  Di(Xo+;AX,) 
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Dj(xo,0 
LDj{xo  +  Axj) 


-Oj(xo) 
£)j(xo  +  Axa) 


U  D](xo  +  Axj ) 


iDj(xo+0  -l)Ax2)  £)2(xo  +  0  -  l)Axj)  6'Z7i(xo  +  (;  -  1  )Ax2) 

[  -Dj(X()  +  jAxj)  Dj{xo+jAxj) 

With  the  further  assumption  that  /,(x)  only  depends  on  the  i"*  component 
of  the  state  vector  we  can  rewrite  (  2.4)  cind  (  2.5)  as 


litJ 

U(x,  t)  =  n  expAi^‘  {E(lAt)  0  U(x,  lAt)]  exp^2A‘  (2.6) 

1=0 

and 


U(x,  t)  =  fl  E(lAt)  0  exp^i^‘  U{x, lAt) exp^2^‘  (2.7) 

1=0 

respectively.  The  symbol  0  is  pointwise  matrix  multiplication  and  I  is  either  an 
j  X  j  or  a.  k  y.  k  identity  matrix.  E{t)  is  defined  as 


o(xo)  a(xo  +  Axi) 

q(x«  +  Axi )  o(xo  +  Axi  +  Axj ) 


o(xo+jAxi)  a(xo -f  jAX) -f  Ax}) 


a(xo  +  tAxj ) 
o(xo  +  Axi  +  tAxj ) 

q(xo  +>Axi  +  r-Axj) 


where  q(x,  t)  =  exp  J— j|l/i||^At  +  C(x)At| .  Note  that  in  discretizing  the  equa¬ 


tions  in  space  the  operation  A"  for  t  =  1  became  matrix  multiplication  on  the 
left  and  for  i  =  2  matrix  multiplication  on  the  right.  If  /,(x)  is  not  assumed  to 
depend  only  on  the  i‘'*component  of  the  state  vector  then  (  2.4)  and  (  2.5)  can 
still  be  implemented  using  matrices,  however,  Aj  changes  for  each  column  of 
the  discretized  density  and  AJ  changes  for  each  row  of  the  discretized  density. 

We  know 


exp^i^^ 


nsO 


(AjAt)^ 

n! 


We  are  interested  in  approximating  the  exponential  of  a  matrix  by  a  truncated 
sum.  If  the  elements  of  AjAt  <  1,  which  we  can  guarantee  by  choosing  At 
small  enough,  then 

exp^i^*  «  I  +  AjAt  . 
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By  using  only  first  and  second  order  approximations  for  the  first  and  second 
order  derivatives  of  the  density  with  respect  to  the  state  we  have  the  row  of  A; 
summing  to  zero.  Additionally,  choosing  At  such  that  the  elements  of  A-,  At  •C  1 
we  have  ensured  the  matrix  I  +  A;At  is  well  conditioned  and  invertible.  We  can 
implement  several  different  schemes  for  (  2.6)  and  (  2.7).  We  can  rewrite  one 
step  of  (  2.6)  as 

U(x,  t  +  At)  =  (1  +  AtA\)  [E(t)  0  U(x,  t)]  (I  +  AtAj )  (2.8) 

which  we  will  denote  as  the  explicit  discretization.  We  can  implement  a  two 
step  approximation  of  (  2.6)  by 

(I- AtA’i)U(x,t  + At)  =  [E(t)0U(x,t)l(H- AtAj)  (2.9) 

U(x,t  +  2At)(I- AtA5)  =  (I  +  AtAJ)(E(t)0U(x,t)|  (2.10) 

which  we  will  denote  as  the  mixed  scheme.  Finally,  we  can  implement  the 
approximation 

(I- AtA*i)U(x,t  + At)(I- AtAz)  =  (E(t)0U(x,t))  (2.11) 

which  we  denote  as  the  implicit  scheme.  (  2.8)  is  called  explicit  because  the 
density  U(x,  t+ At)  is  explicitly  defined  interms  of  the  density  U(x,  t).  (  2.11)  is 
called  implicit  because  each  point  in  the  density  U(x,  t+ At)  is  defined  implicitly 
by  one  or  more  points  in  the  density  U(x,  t  +  At)  as  well  as  U(x,  t)  and  a  set 
of  simultaneous  equations  must  be  solved  to  obtain  U(x,  t  +  At).  The  scheme 
corresponding  to  (  2.9)  and  (  2.10)  having  components  of  both  the  explicit  and 
implicit  schemes  is  denoted  as  the  mixed  scheme.  For  implementing  (  2.7)  we 
again  get  similar  result  with  the  pointwise  multiplication  by  E(t)  done  after 
matrix  multiplication  of  I  +  AtA>  and  U(x,  t). 

Several  approximations  were  used,  in  arriving  at  a  discretization  scheme  for 
(  2.4)  or  (  2.5),  v'hich  warrant  further  discussion.  We  approximated  the  expo¬ 
nential  of  a  matrix  by  using  only  two  terms  of  the  Taylor  expansion.  In  order 
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for  the  approximation  to  be  accurate  it  is  necessary  that  the  elements  of  the 
matrix  be  small.  The  only  term  which  is  free  to  be  set  is  At  but  maliing  At 
small  requires  more  computations  to  process  a  fixed  time  block  of  data.  Al¬ 
ternately,  a  better  approximation  for  the  exponential  can  be  used  by  including 
more  terms  of  the  expansion  but  for  each  additional  term  the  matrices  have 
an  additional  nonzero  upper  and  lower  diagonal  which  again  increases  computa¬ 
tional  complexity  and  Aj  can  nolonger  be  guaranteed  to  be  a  well  conditioned 
matrix.  Another  approximation  used  was  the  first  and  second  order  approxi¬ 
mations  for  the  first  and  second  order  derivatives  of  the  density  u(x,t).  Higher 
order  approximations  can  be  used  but  result  in  additional  nonzero  upper  and 
lower  diagonals  in  AJ  which  again  increases  computational  complexity  and  the 
matrix  I  -I-  AjAt  can  not  be  guaranteed  to  be  well  conditioned.  There  is  also 
a  question  of  which  scheme  to  implement  to  approximate  the  conditional  den¬ 
sity  it(x,t).  The  implicit  schemes  require  more  computations  than  the  explicit 
schemes  but  the  implicit  schemes  exhibit  better  numerical  properties  as  At  be¬ 
comes  large.  However,  since  we  require  At  to  be  small  the  implicit,  explicit  and 
mixed  schemes  arc  essentially  equivalent  with  respect  to  numerical  stability. 
There  is  also  the  approximation  of  the  conditional  density  «(x,  t)by  the  matrix 
U(x,  t).  If  the  discretization  is  very  coarse  U(x,  t)  will  be  a  poor  approximation 
of  the  density  u(x,t)but  the  more  fine  it  becomes  the  more  calculations  per  step 
in  time  are  needed.  Finally,  there  is  a  question  as  to  which  aproximation,  (  2.6) 
or  (  2.7),  is  better  to  implement.  Some  of  these  considerations  have  been  stud¬ 
ied  for  particular  cases  of  interest  to  the  Navy.  These  cases  impose  additional 
constraints  on  the  implementation  scheme  used  and  are  discussed  in  Chapter  3 
and  Chapter  4. 
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2.4  Implementation 


The  discretized  schemes  were  implemented  using  MACSYMA  and  FORTRAN. 
MACSYMA  is  a  symbolic  math  manipulation  paclcage  written  by  the  MIT  Lab¬ 
oratory  for  Computer  Science.  It  is  able  to  generate  legal  FORTR.A.N  code 
which  was  used  to  implement  the  actual  numerical  calculations.  The  use  of 
MACSYMA  gave  significant  flexibility  in  entering  models  by  allowing  symbolic 
calculations  including  differentiation  required  to  implement  the  numerical  solu¬ 
tion.  It  also  allowed  for  easy  changes  to  the  discretization  of  the  state  space. 
After  specifying  the  models  and  state  space  discretization,  FORTRAN  code  is 
automatically  generated  which  simulates  the  diffusion  (using  first  order  differ¬ 
ence  equations)  and  solves  the  Zakai  equation.  The  code  makes  use  of  routines 
from  LINP.A.CK,  a  set  of  FORTRAN  subroutines  for  doing  linear  algebra.  The 
FORTRAN  allows  entry  of  any  constant  parameter  at  the  time  of  execution. 

Presently  the  MACSYMA  code  is  running  on  a  Texas  Instruments  Explorer 
and  the  FORTRAN  code  is  then  transferred  to  a  VMS  or  Unix  machine  and 
executed. 


2.5  Comparison  of  discretization  schemes 

We  presented  five  possible  discretization  schemes  for  the  solution  to  the  Zakai 
equation.  For  (  2.4)  there  are  the  implicit,  mixed,  and  explicit  schemes.  For 
(  2.5)  there  are  the  mixed,  and  explicit  schemes.  Note  that  the  implicit  scheme 
for  (  2.5)  would  be  the  same  as  for  (  2.4).  This  section  gives  results  comparing 
the  numerical  approximation  of  the  Zakai  equation  by  (  2.6)  or  (  2.7)  and  the  use 
of  the  implicit,  mixed  of  explicit  scheme.  A  comparison  for  one  particular  case 
is  shown  in  figures  2.1  to  2.8.  These  figures  correspond  to  a  lognormal  model 
versus  a  Rayleigh  model  with  the  Rayleigh  decorrelation  time  one  tenth  that 
of  the  lognormal,  (reference  tables  4.1  and  4.2).  The  use  of  these  particular 
models  are  discussed  in  Chapter  3.  In  figures  2.1  to  2.4  we  show  the  log 
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likeliliood  ratios  given  hypothesis  0  data  and  in  figures  2.5  to  2.8  we  show  the 
log  likelihood  ratio  given  hypothesis  1  data  for  one  particular  set  of  parameters. 
We  see  that  in  all  cases  the  log  likelihood  ratios  are  very  similar,  .‘^gain,  this 
is  parth’  clue  to  the  small  time  discretization  steps  which  was  necessar}'  for  the 
appro.ximation  to  the  exponential  of  a  matrix  to  be  accurate.  .A.11  the  plots  of  the 
log  likelihood  ratio  are  compared  with  the  implicit  discretization  of  {  2.4)  which 
is  given  by  (  2.11).  From  Monte  Carlo  simulations  presented  in  Chapter  5,  using 
the  implicit  discretization  and  the  mixed  discretization  of  (  2.4),  we  conclude 
that  the  results  are  nearly  identical  in  the  number  of  correct  decisions  and  the 
decision  times  which  supports  the  conclusion  suggested  by  figures  2.1  to  2.8 
that  for  the  given  models  the  approximation  scheme  behave  nearly  identically. 
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log  likelihood  ratio  log  likelihood  ratio 


0 


Log  Likelihood  Ratio  Given  Hypothesis  0 


Figure  2.1:  Mixed  discretization  of  2.4  under  hypothesis  0,  SNR=:5 


Figure  2.2:  Explicit  discretization  of  2.4  under  hypothesis  0,  SNR=5 


log  likelihood  ratio  log  likelihood  ratio 


Figure  2.3:  Mixed  discretization  of  2.5  under  hypothesis  0,  SNR=5 


Figure  2.4:  Explicit  discretization  of  2.5  under  hypothesis  0,  SNR=5 


log  likelihood  ratio  '°8  likelihood  ratio 


Figure  2.5:  Mixed  discretization  of  2.4  under  hypothesis  1,  SNR=5 


Figure  2.6:  Explicit  discretization  of  2.4  under  hypothesis  1,  SNR=5 
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distribution. 

The  typical  radar  system  is  composed  of  an  IF  amplifier  followed  by  a  second 
detector  and  a  video  amplifier.  The  second  detector  and  the  video  amplifier  form 
an  envelope  detector.  If  the  inphase  and  quaxlrature  components  of  the  signal 
entering  the  second  detector  are  independent  and  have  a  Gaussian  distribution 
then  Rice  has  shown  the  density  of  the  envelope  to  be  Rayleigh  (amplitude)  and 
the  power  will  be  distributed  exponentially  (Rayleigh  Power). 

The  Rayleigh  distribution  has  also  been  shown  to  model  sea  clutter  when  the 
resolution  cell  (ie.  the  area  illuminated  by  the  radar  pulse)  is  large  relative  to 
the  water  wavelength.  It  has  also  been  seen  that  sea  clutter  with  high  resolution 
radar  will  often  ha\'e  heavier  tails  which  can  better  be  represented  by  a  lognormal 
density  function. 

In  addition  to  the  amplitude  distribution  of  the  chaff  and  clutter  returns, 
the  power  spectrum  is  also  very  important.  It  gives  information  as  to  the  rate 
at  which  the  radar  cross  section  changes.  Early  models  of  the  spectrum  were 
assumed  to  be  Gaussian  whicli  gives  an  autocovariance  that  is  also  Gaussian. 
More  recent  measurements  performed  in  the  1960’s  and  1970’s  showed  power 
spectrums  of  the  form 

A 


to  be  more  appropriate,  where 

.4  is  the  mean  value  of  the  power  density 
fc  is  the  clutter  spectrum  half-power  frequency 
/  is  frequency 
n  is  a  positive  real  number. 

.According  to  Barton,  the  actual  spectrum  can  generally  be  fitted  with  a  “band- 
limited  white  noise  spectrum,  extending  from  zero  to  infinity,  but  with  steadily 
decreasing  amplitude  above  (a)  given  frequency”  (Bar64,p.82|,  which  we  denote 
by  fe  (the  half-power  frequency).  Then  the  actual  spectrum  can  be  well  approx- 
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imated  by 


I-''  (a) 

The  autocovariance  function  is  given  by 

R{t)  =  /l7r/cexp{-27r/c|tl} 

which  results  in  exponential  decorrelation  with  time  constant  5^. 

In  the  case  of  a  ship  the  scattering  elements  are  much  more  complex  and  most 
of  the  modelling  has  been  based  on  empirical  considerations.  “The  lognormal 
distribution  has  been  used  to  model  scattering  from  highly  directive  reflectors 
when  viewed  from  random  aspects”  [Bla86,p.77],  such  as  randomly  oriented  flat 
plates,  corner  reflectors,  and  antennas.  “It  has  become  popular  to  use  the  lognor¬ 
mal  model  to  describe  large  metal  objects  of  irregular  shape,  especially  ocean 
vessels  with  complex  superstructures”  [MM73,p.67].  We  do  not  wish  to  imply, 
however,  that  the  lognormal  distribution  can  model  all  RCS  returns  from  ships. 
In  addition  to  the  lognormal  distribution,  radar  cross  section  (RCS)  returns  have 
also  been  observed  to  be  distributed  Rayleigh  power,  Rayleigh  amplitude,  chi- 
square  with  two  and  four  degrees  of  freedom,  and  Rician.  Lognormal  statistics 
tend  to  appear  at  major  aspects  of  ships  with  dominant  scatterers,  especially 
on  the  larger  ships  such  as  carriers  and  battleships. 

As  in  the  case  of  chaff  and  clutter,  it  has  been  observed  “that  the  power 
spectral  density  of  the  ship  RCS  fluctuation  due  to  deterministic  and  random 
azimuth(aspect),  pitch,  and  roll  motion  can  be  very  well  approximated  by  a 
power  spectral  density  of  the  form”  [Bar80,p.30] 

which  yields  the  autocovariance  function 

R(<)  =  A7r/iexp{-27r/tltJ} 

which  is  exponential  decorrelation  with  time  constant  2^-  It  has  been  observed 
that  the  decorrelation  time  for  a  ship  tends  to  be  larger  than  for  chaff  or  clutter 
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Chapter  3 

Model  Definition  and 
Parameter  Estimation 


3.1  Specification  of  model 

We  are  interested  in  the  performance  of  the  detector  for  the  previously  described 
discretization  schemes  and  as  such  we  need  specific  models.  The  Navy  is  inter¬ 
ested  in  the  application  of  the  resulting  detector  to  radar  problems.  Based  on 
these  radar  problems  we  choose  specific  signal  and  noise  models.  We  will  look 
at  the  ship  versus  decoy  problem. 

One  common  decoy  used  is  chaff.  Chaff  is  the  code  name  used  during  WW  II 
to  refer  to  metallic  dipoles  (strips  of  light  metal  foil)  dropped  or  launched  into  the 
air  to  confuse  enemy  radar  by  presenting  a  large  radar  return.  If  these  dipoles 
are  half  the  length  of  the  radar  wavelength  they  will  resonate  and  give  large 
radar  returns.  Since  the  chaff  cloud  consists  of  a  large  number  of  scatterei-s 
of  which  no  single  one  dominates,  the  phase  and  amplitude  variations  of  the 
individual  scattering  elements  can  be  considered  independent.  By  application  of 
the  Central  Limit  theorem  to  the  phasor  summation  of  the  RF  voltages  induced 
in  the  receiving  antenna  by  the  reflection  of  the  individual  scatterere  it  leads  to 
the  conclusion  that  the  inphase  and  quadrature  components  have  a  Gaussian 
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which  is  supported  by  the  fact  that  the  individual  scattering  elements  of  a  ship 
are  fixed  relative  to  each  other  and  the  aizimuth,  pitch,  and  roll  motion  is  slow 
relative  the  orientation  changes  of  the  dipoles  in  the  falling  chaff  or  the  sea 
surface. 

The  Rayleigh  and  lognormal  radar  returns  can  be  modeled  using  diffusions 
which  are  given  by 

dxt  =  r{xt)dt  +  g'{xt)dW\ 
dyt  =  h'{xt)dt  + 'rdV t  • 

We  will  drop  the  time  subscript  and  use  subscripts  to  denote  elements  in  vectors 
or  matrices.  Under  hypothesis  0  we  will  assume  that 

/®(x)  =  a 

b  O' 

0  b 

h°{x)  =  (c(xf +  x|)*], 

then  /i°(x)  has  Rayleigh  amplitude  statistics.  As  for  hypothesis  1,  we  assume 
that 

/‘(x)  =  [r/xi  ] 

=  [r] 

/i‘(x)  =  (sexp{xi}] , 

then  h^{x)  is  distributed  according  to  a  lognormal  density.  The  parameters  a, 
6,  c,  q  ,r,  s,  and  7  are  constants.  Tlie  term  7cfV ,  models  thermal  noise  in  the 
receiver.  The  constant  7  must  be  nonzero  or  else  the  likelihood  ratio  A,  becomes 
undefined  which  results  in  numerical  problems  in  the  discretization. 

We  will  wish  to  test  the  detector  with  known  parameter  \'alues  and  with 
estimated  values.  The  parameters  we  need  to  estimate  are  the  constants  a,  i,  c, 
7  in  the  Rayleigh  case  and  q,  r,  s,  7  in  the  lognormal  case.  In  the  case  7 


Xi 

X2 
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as  explained  earlier,  the  term  7  becomes  a  renormalizing  term.  When  a  and  b 
[q  and  j  )  are  independent  of  time  the  underlying  state  is  a  correlated  Gaussian 
path.  For  the  discretized  version  of  the  differential  equation  the  parameters  a 
and  b  (q  and  r)  can  then  be  related  to  the  parameters  pc  the  correlation  and 
(Jg  stationary  variance  of  the  underlying  Gaussian  (appendix  D)  by 


'■“"H — s — ) 

where  pc  =  exp  and  to  is  the  decorrelation  time  constant  for  the  under¬ 

lying  Gaussian. 

3.2  Parameter  estimation:  lognormal 

One  very  common  estimator  is  the  maximum  likelihood  estimator  (MLE).  To 
calculate  the  MLE  of  the  parameters  for  the  discretized  version  of  the  solution 
to  the  Zakai  equation  we  need  to  find  the  distribution  of  the  output  data.  To 
generate  the  lognormal  or  Rayleigh  distribution  given  the  form  of  the  models  we 
firet  look  at  the  distribution  for  the  underlying  Gaussian.  With  the  requirement 
that  ^  <  a  <  0  (or  q  respectively),  where  At  is  the  time  discretization  used, 
then  the  underlying  Gaussian  has  a  stationary  distribution  with  zero  mean, 
variance  denoted  cr^-,  and  exponential  decorrelation  ic.  R(k)=  exp  •c  =  p^. 
We  can  write  the  distribution  for  the  underlying  Gaussian  in  closed  form  for 
both  the  lognormal  and  the  Rayleigh  cases,  refer  to  appendix  A. 

For  the  lognormal  case  /i*(x)  is  an  invertible  transformation,  so  we  can  use 
the  Jacobian  change  of  variables  formula  to  determine  the  distribution  for  /i‘  (x). 
Denote  y,  =  /i*(x,)  then  we  can  write  the  distibution  for  y  =  [yi,y2,  ••■,yn]  as 
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where 


exp 


n— 1  n-l  y 

0  rtj^  \  '  fil  fi  -4-  ll  4”  \  ^  f/1” 

“Pc  ij  Pc 

,=i  1=2  JJ 

[i]  =  In  ^  -  n  =  In(yi)  -  (ln(s)  +  ^i)  . 

To  calculate  the  distribution  of  the  observation  (h'(x)  plus  additive  noise) 
in  closed  form  is  not  possible  so  one  cannot  find  the  MLE,  however,  if  the  signal 
to  noise  ration  (SNR)  is  very  large  the  ML  estimator  derived  assuming  no  noise 
may  be  a  good  estimator.  The  ML  estimators  , assuming  no  noise, for  <75  and 
ln(s)  —  //  with  respect  to  pc  are  given  by 


ln(s)  -  /t  = 


In(tyi)  +  In(yn)  +  (1  -  Pc)  in(7/,) 
2  +  {n-2)(l-pc) 


where 


= 


«(1  -  Ph)  L 


n— I 


n  —  l 


-  2pG  +  11  +  Pc  E  W" 


1=1 


t=l 


1=2 


[i]  =  In(yi)  -  (ln(s)  +  p)  . 


The  ML  estimator  for  pc  cannot  be  solved  for  explicitly,  however,  we  are 
only  interested  in  solutions  in  the  interv'al  (—1,1)  and  can  easily  implement  an 
algorithm  to  find  the  roots  of  the  derivative  of  the  density  with  respect  to  pc 
in  this  interval.  Differentiating  the  lognormal  density  with  respect  to  pc  and 
setting  it  equal  to  zero  we  get 


0  =  (ir  -  l)<T^Pc(l  ~  Pc)  ~  Pg  [^1  —  2pcT2  +  PcTaj  +  (1  —  Pc)  [^2  ~  P^s]  (3-1) 


with 


Tx 

Ti 

Tz 


i=i 


XI  ln-(i/.)  -  2p  ^  ln(y<)  +  ny? 

£bi(T/.)  +  flln(?/,) 

=1 


1=1 

n— 1 

Y^lii(t/,)ln(7/i+i)  -Ji 

tsl 

£  ln-(T/<)  -  2p  ln(y,)  +  (n  -  2)p' 

1=2  is2 


1=2 

772 


(n  -  l)p' 
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which  depend  on  pg  thru  JI,  the  ML  estimate  for  ln(5)  +  A*  for  ^  particular  value 
of  Pc-  D\'  finding  the  root  of  (  3.1)  we  obtain  an  estimate  for  pc  and  in  turn  an 
estimate  for  ln(s)  +  p  and  Cc-  Note  that  the  ML  estimator  for  can  also  be 
written  as 

"  n(i-;4)  • 

These  estimates  are  not  very  good  when  any  significant  amount  of  additive 
white  noise  in  the  output  equation  is  present,  particularly  for  pc-  As  will  be 
seen  later  the  detector  becomes  numerically  unstable  if  the  SNR  is  too  large, 
ie.  small  noise  power  relative  to  signal  power,  so  the  ML  estimates  assuming  no 
noise  are  not  useful  in  this  case.  An  alternate  estimation  scheme  would  be  to 
estimate  pc  and  then  use  the  ML  estimator  for  correlated  lognormal  data  or  the 
ML  estimators  for  independent  data  which  aie 

ln{s)  +  p= 


al  = -  {\n{s)  +  p)  - 

”  i=i 

We  can  estimate  pc  from  the  covariance  function  of  the  output  data.  If  the 
underlying  state  is  Gaussian  with  an  exponential  covariance  function  then  the 
lognormal  sample  path  hits  a  covariance  function  which  is  approximately  expo¬ 
nential  (DarSO).  Assuming  an  exponential  covariance  function  for  the  lognormal 
signal  we  match  an  exponential  to  Ro{k)  ,k  >  0(o  denoting  observation  data). 
From  this  we  get  an  estimate  for  the  decorrelation  time  of  the  lognormal  data 
which  we  will  denote  as  t(.  Then  we  use  the  approximation 


ln(a5)  -  In 


In 


exp(a^)  -l\Y 

exp(l)  /J. 


from  [BarSO]  to  get  an  estimate  for  tc  and  thus  pc-  We  justify  the  estimate 
for  ti  by  the  fact  that  the  lognormal  data  without  additive  white  noise  has  a 
covariance  function  which  is  approximately  exponential ,  the  white  noise  covari¬ 
ance  function  is  a  Delta  function,  and  since  the  lognormal  data  and  the  white 
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noise  are  independent  the  covariance  function  of  the  observation  is  the  sum  of 
the  covariance  function  of  the  lognormal  signal  and  the  co\'ariance  function  of 
the  noise.  Figures  3.1.  3.3  and  3.5  are  plots  of  lognormal  model  sample  paths 
without  noise  and  with  additive  noise  for  SNR  of  20,  5  and  2.  Figures  3.2,  3.4 
and  3.C  are  plots  of  the  normalized  covariance  functions  of  the  signal  without 
noise  (solid),  the  signal  plus  noise  (short  and  long  dash)  Ro{k),k  >  0,  signal 
plus  noise  with  the  first  term  removed(short  dash)  Roik),  k  >  0,  and  an  e.\po- 
nential  function  (long  dash).  In  figures  3.10  to  3.15  we  have  similar  plots  for 
the  Rayleigh  model  for  SNR=5  and  decorrelation  times  U  =  see 

that  the  covariance  function  does  exhibit  the  desired  behavior  of  a  6  function 
at  zero  for  the  additive  noise  plus  an  approximately  exponential  function  for 
Ro{k),k  >  0.  .\n  additional  observation  is  that  as  the  SNR  becomes  smaller 
(more  noise)  the  exponential  behavior  of  the  covariance  function,  Ro{k),k  >  0, 
becomes  more  corrupted  due  to  the  noise  covariance  function  not  being  a  true 
6  function,  esnecially  for  the  lognormal  case  where  the  stationarj'  variance  is 
smaller  than  for  the  Rayleigh  case  and  the  signal  begins  to  look  like  a  mean 
which  is  removed  by  the  covariance  function. 

Another  possible  estimator  for  ln(s)+;i  and  would  be  the  mean  to  median 
ratio  of  i.i.d.  data  (reference  appendix  B).  If  we  assume  the  underlying  state 
has  zero  mean  then 

mean 

— - —  =  exp 
median  2 

median  =  s  exp  ^  =  s  . 

Finally,  from  Huber  we  have  a  robust  estimator  for  independent  identically 
distributed  Gaussian  data  known  as  the  Median  ,'\bsolute  Deviation  (MAD) 
estimator  [HubSlj.  By  taking  the  logarithm  of  independent  noisy  data  values 
(rejecting  any  negative  obseiA'ations)  we  have  corrupted  iid  Gaussian  data.  Let 
r,  =  In  y,  then 

ln(5)  +  ft  =  median{z,} 

Oq  =  [median{lz,  -  median{2,}|}]* 
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We  again  estimate  pc  from  the  covariance  function. 


3.3  Parameter  estimation:  Rayleigh 


For  the  Rayleigh  ciise  if  we  set  y,  =  h°(x,)  we  have  the  distribution  of  y  = 
as 


/(y)  =  ^2- (n-l) 


II  ^olPyiVi+i] 

i=l 


Ltsl 


where 


1  =  2  71  —  1 

■leJn-o-L>  i  —  1 


a,  =  ^  2<.fcc7a-Pc) 


i  =  \.n 


PG 


^  a};c\l-p-G) 

and  Iq  is  the  modified  zero  order  Bessel  function.  It  is  not  possible  to  estimate 
c  and  og  independently  since  c  and  oq  alway  appear  together.  We  can  only 
estimate  the  product  cog.  Since  we  can  only  estimate  the  product  of  c  and  oq 
we  arbitrarily  set  the  stationary  variance  of  the  Rayleigh  model  to  the  stationary 
variance  of  the  lognormal  model  when  we  implement  the  detector. 

The  ML  estimator  for  cac  and  pc  is  the  solution  to  the  equations 

0  =  -inc'c-Gi  1  -  P5)  -  +  yi;  +  { 1  +  /?■)  £  i/i') 


n— 1  Ij 

•^PGH^yiVi+v  ( 


and 


0  =  2pG{n  -  1 )( 1  -  -  (y'l  +  y;  +  2  ^  yf  j 

Ii  (-^ff^y.y.^i) 

+(i  +/?g)  - r 

•=*  lo  (?.gfit4)y.y.+ij 

where  /j  is  the  modified  first  order  Bessel  function.  Because  of  the  poor  results 
for  the  ML  estimator  for  the  lognormal  case  when  noise  is  presen .  and  the  fact 
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Figure  3.3:  Lognormal  sample  path,  SNR=5 


Figure  3.4:  Covariance  functions  for  Lognormal  sample  path,  SNR=5 


Figure  3.6;  Covariance  functions  for  Lognormal  sample  path,  SNR=2 


that  the  detector  becomes  numerically  unstable  for  large  SNR  (ie.  small  additive 
noise  power)  this  estimator  was  not  implemented. 

An  alternate  method  of  estimating  the  Rayleigh  parameters  is  to  estimate 
Pc  from  the  covariance  function  of  the  output  data  similar  to  what  we  did  in 
the  lognormal  case.  If  xji  is  correlated  Rayleigh  (amplitude)  data  and  we  define 
-i  =  y]  then  r,  is  correlated  Rayleigh  power  data.  If  the  underlying  state  is 
Gaussian  with  exponential  decorrelation  and  time  constant  tQ  then  r,  will  have 
exponential  decorrelation  with  time  constant  f.  =  [BarSO].  Figures  3.7  to 
3.9  show  the  normalized  covaricince  functions  of  a  Rayleigh  amplitude  signal,  the 
square  of  the  Rayleigh  amplitude  signal  and  a  matched  exponential  function.  We 
see  that  the  decorrelation  time  tr  «  t..  We  note  that  as  the  decorrelation  time 
increases, Pc,"  1  •  the  covariance  function  does  not  go  smoothly  to  zero.  The 
oscillatory  behavior  is  due  to  numerical  instabilities  in  the  discretization  scheme 
as  Pc  -*  1  and  is  known  as  ringing.  Figures  3.10,  3.12  and  3.14  are  plots  of 
Rayleigh  (amplitude)  model  sample  paths  without  noise  and  with  additive  noise 
for  SNR  5.  Figures  3.11,  3.13  and  3.15  are  plots  of  the  normalized  covariance 
functions  of  the  signal  without  noise,  the  signal  plus  noise  Ro{k),k  >  0,  signal 
plus  noise  with  the  first  term  removed  Ro{k),  k  >  0,  and  an  exponential  function. 
As  noted  earlier  wn  ha\'c  the  desired  behavior  that  the  covariance  function  of 
the  ob.servation  is  approximately  the  sum  of  an  exponential  function,  due  to  the 
signal,  and  a  6  function  due  to  the  additive  noise. 

We  can  estimate  ede  from  the  ML  estimator  for  ctcq  for  i.i.d.  observations 
which  is  given  by 


Alternately,  we  could  use  the 


2  2  ^  *> 


power  to  median  or  the  power  to  mean  ratios  for 


i.i.d.  obserr'ations  given  by  (reference  appendix  C) 


power  _  2 

median  “ 
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power 

mean 


8 

-cac 
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3.4  Estimation  of  noise  scaling  7 

We  will  estimate  7  from  the  covariance  function  of  the  observations.  The  noise 
in  the  output  equation  is  white  zero  mean  Gaussian  so  its  covariance  function 
is  a  6  function.  Since  the  noise  is  independent  of  the  signal  /i(x)  the  covariance 
function  of  the  oijservation  data  is  the  sum  of  the  covariance  functions  of  the 
noise  and  the  signal.  We  use  R{k)  to  denote  the  covariance  function,  with 
subscripts  s.n,o  to  denote  signal,  noise,  and  observation  respectively.  Under 
the  model  assumptions  made  the  signal  covariance  function  is  approximately 
exponential  so  71,(0)  is  approximately  7?,(1)  provided  the  decorrelation  time 
constant  of  the  signal  is  not  on  the  order  of  the  time  discretization  At.  With 
the  noise  covariance  function  assumed  to  be  a  ^  72,(1)  «  72o(l)-  So  we  have 

72„(0)  72,(0)  +  72n(0) 

or 

72„(0)  =  72„(0)  -  72,(0)  . 

So  72,(0)  w  72„(1)  and  substituting  we  get 

72„(0)5s72„(0)-72o(1)  . 

In  figures  3.11  to  3.15  as  the  decorrelation  time  approaches  the  time  discretiza¬ 
tion  step  size  At  the  6  function,  due  to  the  noise,  becomes  less  distinct  because 
the  covariance  function  of  the  signal  is  approaching  a  6  function  with  respect 
to  the  time  step  size  At.  Additionally,  if  the  signal  to  noise  ratio  is  large,  ic. 
the  noise  covariance  function  has  a  6  function  with  magnitude  on  the  order  of 
72,(0) -72,(1),  then  72„(0)  «  72o(0)  — 72„(1)  is  clearly  a  poor  approximation  which 
could  cause  significant  numerical  problems  and  errors  since  the  observations  are 
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Figure  3.7:  Covariance  functions  for  Rayleigh  amplitude  and  Rayleigh  power 
sample  path,  tr  100 
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Figure  3.S:  Covariance  functions  for  Rayleigh  amplitude  and  Rayleigh  power 
sample  path,  tr  ~  t//10 
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Figure  3.9:  Covariance  functions  for  Rayleigh  amplitude  and  Rayleigh  power 
sample  path,  <r  ^  h 
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Figure  3.10:  Rayleigh  sample  path,f,  «  t//100  SNR=5 


Figure  3.11:  Covariance  functions  for  Rayleigh  sample  path,  SNR=5 
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renormalized  by  ^  where  i2„(0)  =  ^  after  discretizing.  Again,  for  the  models 
we  use.  figures  3.11,  3.13  and  3.15  exhibit  the  desired  behavior  of  a  (5  function 
plus  an  exponential  function.  This  estimate  could  be  improved  by  extrapolating 
back  to  /2o(0)  after  matching  the  exponential  to  Ro{k)  k  >  0.  This  would  help 
remove  variations  of  Roik)  from  Rs(k)  for  A:  >  0  due  to  the  additive  noise  not 
being  an  exact  6  function  as  seen  most  significantly  in  figures  3.2,  3.4,and  3.6. 
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Chapter  4 


Results 

4.1  Model  specification 

We  choose  the  Rayleigh  and  lognormal  models  based  on  NRL  data.  Tables  4.1 
and  4.2  summarize  the  parameter  values  chosen  to  test  the  detectors  perfor¬ 
mance.  Results  from  Chapter  3  and  the  appendices  were  used  to  determine 
parameter  values  and  signal  and  noise  statistics.  The  Rayleigh  model  with  pa¬ 
rameters  A  =  -24.06, B  =  1.038  and  C  =  .2117  corresponds  to  the  Chaff/decoy. 
We  will  denote  the  hypothesis  test  between  this  particular  Rayleigh  model  and 
the  lognormal  model  as  the  ship  versus  decoy  case.  The  other  se?^f  parame- 
tere  are  to  test  the  effect  of  different  characteristics  of  the  signal  and  noise  on 
the  detector  performance.  The  decorrelation  time  and  the  stationary  mean  are 
investigated.  The  second  set  of  five  parameters  for  the  Rayleigh  have  the  sta¬ 
tionary  mean  shifted  by  5  percent  and  the  powers  moved  closer.  In  addition  to 


Lognormal  parameters 

Q  R  S 

stationary  statistics 

(Tq  (ii  af  power 

-4.800  .4646  .03935 

.02251  .03980  .00003605  .001620 

.2066 

Table  4.1:  Lognormal  model  parameters 
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Rayleigh  parameters 

ABC 

stationary  statistics 

(Tq  Hr 

power 

-.04840 

.04668 

.2117 

.02251 

.03980 

.0004330 

.002018 

lOOtct 

-.2420 

.1044 

.2117 

.02251 

.03980 

.0004330 

.002018 

lOtci 

-2.419 

.3299 

.2117 

.02251 

.03980 

.0004330 

.002018 

tci 

-24.06 

1.038 

.2117 

.02251 

.03980 

.0004330 

.002018 

tci /lO 

-227.9 

3.111 

.2117 

.02251 

.03980 

.0004330 

.002018 

tcil  100 

-.04840 

.04668 

.2011 

.02251 

.03781 

.0004330 

.001821 

lOOtci 

-.2420 

.1044 

.2011 

.02251 

.03781 

.0004330 

.001821 

lOtci 

-2.419 

.3299 

.2011 

.02251 

.03781 

.0004330 

.001821 

tc, 

-24.06 

1.038 

.2011 

.02251 

.03781 

.0004330 

.001821 

tc{/l0 

-227.9 

3.111 

.2011 

.02251 

.03781 

.0004330 

.001821 

tci/ 100 

Table  4.2:  Rayleigh  model  parameters 

looking  at  different  values  for  a,  b,  c,  q,  r  and  s  we  also  investigate  the  detec¬ 
tors  Ijchavior  for  different  SNR.  Table  4.3  summarizes  the  values  of  7  and  the 
corresponding  SNR.  We  also  tested  the  detectors  behavior  for  different  values 
of  Pp  and  P.\t  though  for  most  simulations  Pp  =  Pm  =  .001  .  We  made  a  more 
indepth  comparison  of  the  implicit  and  the  mixed  schemes  implementing  (  2.4) 
than  it  was  done  in  Chapter  1.  The  probability  of  detection,  the  probability  of 


7 

power 

SNR  w.r.t  stationary 

Lognormal  power 

.0006364 

.000810 

2 

.0004025 

.000324 

5 

.00020125 

.000081 

20 

Table  4.3; 

Noise  scale 
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miss  and  the  average  detection  time  are  calculated  from  1000  simulation  runs. 


4.2  Simulations 

The  following  section  summarizes  results  of  simulations  for  the  models  specified 
in  tables  4.1  and  4.2  for  the  SNR  given  in  table  4.3.  For  these  simulations 
Pp  =  =  .001.  Effects  of  the  decorrelation  time,  SNR,  stationary  mean  and 

implicit  \-ersus  mixed  scheme  on  the  percentage  of  correct  decisions  and  average 
decision  time  are  given.  The  results  are  based  on  1000  runs.  Figures  4.1  to  4.8 
give  results  for  the  percentage  of  correct  decisions  versus  the  decorrelation  time 
(of  the  Rayleigh  relative  to  the  lognormal)  under  hypothesis  0  and  1.  The  plots 
are  for  SNR  of  20,  5  and  2,  with  the  (stationary)  means  matched  and  shifted  by 
5  percent,  and  for  the  implicit  and  mixed  scheme  implementing  (  2.4).  Figures 
4.9  to  4.16  give  the  corresponding  plots  for  the  average  detection  time  versus 
decorrelation  time.  These  plots  correspond  to  tables  E.l  to  E.6  in  appendix  E. 

4.2.1  Effects  of  decorrelation  time,  SNR,  mean,  and 

scheme  on  percentage  of  correct  decisions 

Under  hypothesis  1  (lognormal)  the  detector  performs  well  with  respect  to  tlie 
percentage  of  correct  decisions  except  for  the  case  tr  ~  t//100  .  In  this  case 
the  detector  performs  best  for  SNR=20,  and  progressively  worse  for  SNR=5 
and  2.  One  possible  reason  for  this  behavior  has  already  been  mentioned  in 
chapter  3;  that  is  as  the  SNR  decreases  (more  noise)  the  covariance  function  of 
the  observation  for  the  lognormal  tends  toward  a  6  function  with  the  lognormal 
signal  behaving  like  a  mean  for  the  noise  and  being  removed.  In  the  cases  U  ~ 
f//100  and  f//10  the  theoretical  decorrelation  time  of  the  Rayleigh  is  approaching 
the  time  step  size  At  and  the  covariance  function  of  the  signal  is  becoming  a  6 
function.  So  we  see  that  the  covariance  function  of  the  lognormal  plus  noise  is 
approaching  the  covariance  function  of  the  Rayleigh.  This  behavior  is  shown  in 
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Percentage  Correct  Percentage  Correct 


Decorrelation  Time  (  scale  log(tT/t,)  ) 


Figure  4.5:  Percentage  of  correct  decisions  under  Hypothesis  1 


Decorrelation  Time  (  scale  log(tr/t,)  ) 

Figure  4.6:  Percentage  of  correct  decisions  under  Hypothesis  1 
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Percentage  Correct  Percentage  Correct 


Decorrelation  Time  (  scale  log(4/li)  ) 

Figure  4.7;  Percentage  of  correct  decisions  under  Hypothesis  1 


Decorrelation  Time  (  scale  log(t,/l,)  ) 


Figure  4.8:  Percentage  of  correct  decisions  under  Hypothesis  1 
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figures  3.1  to  3.6  and  figures  3.10  to  3.15. 

Under  lu'pothesis  0.  since  this  is  a  binaiy  hypothesis  test,  as  e.xpected  we  see 
the  detector  works  well  for  U  ss  i//100  and  t//10  and  does  poorly  for  U  ~  t/.lOt/ 
and  loot/.  We  also  have  the  SNR  effecting  the  percentage  of  correct  decisions 
opposite  its  effect  under  hypothesis  1.  This  may  be  due  to  the  fact  that  for  a 
higher  SNR  (less  noise,  smaller  7)  the  observations,  which  are  renormalized  by 
are  weighted  more  in  the  solution  to  the  Zakai  equation  and  cause  a  quicker 
decision  liut  for  these  longer  decorrelation  times  fewer  uncorrelated  blocks  of 
data  are  ob.servcd.  The  smaller  SNR  allows  a  longer  observation. 

Under  both  lypothesis  we  see  that  shifting  the  mean  had  only  a  small  effect. 
Under  hypothesis  1  there  was  a  decrease  in  performance  with  respect  to  the 
percentage  of  correct  decisions  while  under  hypothesis  0  there  was  a  slightly 
more  significant  increase  in  performance.  Several  runs  were  made  when  the 
means  were  significantly  different  (a  factor  of  10)  and  one  could  easily  distinguish 
the  two  hypothesis.  In  this  case  the  detector  made  correct  decisions  very  quickly. 
However,  when  the  means  are  so  different  these  cases  are  not  of  interest  since 
much  simpler  detectors  could  be  implemented.  With  regards  to  the  two  scheme 
implemented,  the  implicit  and  mixed,  for  (  2.4),  we  see  that  the  percentage  of 
correct  decisions  is  nearly  identical. 

4.2.2  Effects  of  decorrelation  time,  SNR,  mean,  and 
scheme  on  average  detection  time 

Under  both  hypothesis  1  and  hypothesis  0  we  see  that  for  small  SNR  (more 
noise)  the  detector  takes  longer  to  make  a  decision  which  one  would  expect. 
This  is  due  to  the  renormalization  of  the  observ’ations  by  As  7  increases 
(smaller  SNR,  more  noise)  the  observations  are  weighted  less  in  the  solution  to 
the  Zakai  equation. 

We  also  see  that  the  mean  has  a  small  effect  on  the  average  detection  time. 
Under  both  hypothesis  the  shifted  mean  usually  increase  detection  time  slightly. 
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Figure  4.9:  Average  decision  time  under  Hypothesis  0 
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Figure  4.10:  Average  decision  time  under  Hypothesis  0 
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Figure  4.11:  Average  decision  time  under  Hypothesis  0 


Figure  4.12:  Average  decision  time  under  Hypothesis  0 
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Figure  4.17:  Average  decision  time  for  different  values  of  Pf  and  Pm 

Again,  the  implicit  and  mixed  scheme  give  nearly  identical  results. 

It  is  the  decorrelation  time  that  has  the  most  significant  effect  on  the  average 
detection  time.  Under  hypothesis  1  for  U  «  t//100  we  had  no  detections  for 
SNR=5  and  2  so  no  average  detection  times  are  plotted.  Except  for  hypothesis 
1  with  SNR=2  we  see  the  detection  time  is  peaked  at  tr  ^  U  and  decreases  as  U 
increases  or  decreases.  This  implies  the  decorrelation  time  is  a  very  significant 
feature  of  the  signal  for  the  detector.  When  the  decorrelation  times  are  matched 
the  detector  takes  the  longest  to  reach  a  decision.  We  also  see  that  for  the 
“symmetric”  cases,  ie.  tr  «  t//100  and  tr  a  lOOt/  or  tr  w  t//10  and  tr  w  lOt;, 
that  the  detection  time  for  the  longer  decorrelation  times  is  longer. 


4.3  Performance  varying  Pp  and  Pm 

In  figure  4.17  we  have  a  plot  of  the  average  detection  time  for  different  values 
of  Pf  and  P^/  for  the  ship  versus  decoy  case  with  the  SNR=5.  The  figure 
corresponds  to  tables  E.7  and  E.8  in  appendix  E.  We  see  that  as  the  Pp  or 
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Lognormal  (Hypothesis  1) 

q 

r 

s 

7 

actual 

-4.80 

.4646 

.03935 

.000403 

estimated 

-6.06 

.4181 

.03933 

.000403 

Rayleigh  (Hypothesis  0) 

a 

b 

c 

7 

actual 

-24.06 

1.038 

.2117 

.000403 

estimated 

-24.06 

2.098 

.1203 

.000410 

Table  4.4:  Estimated  parameters  for  the  ship  versus  decoy  case,  SNR=5 

the  P,\i  increases,  for  fixed  P\{  and  Pp  respectivly,  the  detection  time  decreases. 
This  is  as  one  would  expect  since  increasing  Pp  or  Pm  allows  greater  error.  VVe 
see  in  tables  E.7  and  E.8  that  the  actual  percentage  of  correct  decisions  is 
usually  better  than  the  theoretical  values  used  to  set  the  thresholds.  For  the 
cases  where  Pp  =  .0001  and  P^  =  .0001  the  percentage  of  correct  decisions  is 
not  significant  since  it  is  base  on  only  1000  runs,  however,  the  average  detection 
time  still  is. 


4.4  Performance  with  estimated  parameters 

In  chapter  3  we  gave  some  justification  for  the  models  we  chose  to  study  the 
detector,  as  well  as  some  techniques  for  estimating  the  parameters  for  these 
models.  We  applied  some  of  these  techniques  to  the  ship  versus  decoy  case 
and  ran  the  detector  with  estimated  parameters  to  determine  what  degradation 
would  occur.  For  a  better  comparison  of  the  results,  the  sample  paths  of  the 
runs  were  the  same  as  those  used  for  the  ship  versus  decoy  case  with  SNR=5 
and  implicit  scheme. 

To  estimate  the  parameters  for  the  lognormal  case  we  generated  a  correlated 
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Hypothesis  0 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  dot. 

time  wrong 

avg.  dot. 

time 

known 

100.0 

68.2 

0 

- 

68.2 

estimated 

100.0 

75.8 

0 

- 

75.8 

Hypothesis  1 

%  coiTect 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

known 

99.5 

393.0 

0.5 

164.4 

391.9 

estimated 

99.6 

294.1 

0.4 

159.5 

293.5 

Table  4.5:  Comparison  of  ship  vei-sus  decoy  case  with  known  and  estimated 
parameters  for  the  implicit  scheme,  SNR=5 


sample  path  with  SNR=5.  We  match  an  exponential  function  to  the  covari¬ 
ance  function  of  the  obsen'ation,  Ro(^)  >  0,  and  estimated  the  lognormal 

decorrelation  time  constant  to  be  .16  seconds.  We  used  the  approximation 


to 


In(ac)  -  In 


In  1-1- 


exp((7^)  -  1 
e.xp{l) 


to  estimate  the  decorrelation  time  of  the  underlying  Gtiussian.  We  related  tlie 
decorrelation  time  constant  of  the  underlying  gaussian  to  the  diffusion  parameter 
q  using  results  from  Appendix  D  to  arrive  at  the  estimate  q  =  —6.06.  From 
independent  identically  distributed  data  we  used  the  median  to  estimate  the 
scale  parameter  s  to  be  .03933  .  Using  the  mean  to  median  ratio  we  estimated 

to  lie  .01444  w’hcih  we  related  to  the  parameter  r,  using  results  in  appendix 
D,  to  get  the  estimate  .4181  for  r. 

To  estimate  the  Rayleigh  parameters  we  again  generate  a  correlated  sample 
path  with  SNR=5.  We  matched  an  exponential  function  to  the  covariance  func¬ 
tion  of  the  obser\'ation,  Ro{k)  k  >  0,  and  estimated  the  Rayleigh  decorrelation 
time  constant  to  be  .019  seconds.  So  the  underlying  Gaussian  has  an  approx- 
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imate  decorrelation  time  constant  of  .038  seconds.  Using  appendix  D  we  got 
an  estimate  of  -24.056  for  a.  From  independent  identically  distributed  data  we 
used  the  a\'erage  of  the  power  to  mean  ratio  which  yielded  an  estimated  cerc  of 
.03773  and  the  power  to  mediam  ratio  which  resulted  in  an  estimate  of  .03645 
to  arrive  at  the  estimate  cffa  =  03709.  As  noted  in  chapter  3  since  c  and  aa 
always  appear  together  in  the  density  we  can  only  estimate  there  product  so  we 
arbitrarily  set  for  the  Rayleigh  model  equal  to  that  of  the  lognormal  model. 
Thus  using  appendix  D  we  estimated  b  to  be  2.098  and  the  scale  parameter  c 
to  be  .1202.  To  get  a  better  idea  of  the  performance  of  the  estimates  for  the 
Rayleigh  model  if  we  set  the  scale  parameter  c  to  its  actual  value  of  .2117  then 
the  estimated  value  of  6  becomes  1.2122  which  appears  to  be  a  much  better  esti¬ 
mate.  Clearly,  since  ccc  appears  together  in  the  correlated  Rayleigh  density  the 
signal  statistics  are  unchanged  regardless  of  how  one  assigns  the  values  of  c  and 
<tg  provided  their  product  remains  the  same.  Clearly  the  conditional  density 
ii(x,t)which  is  the  solution  to  the  Zakai  equation  would  be  different  though  the 
statistics  of  the  signal  are  theoretically  the  same.  However,  due  to  the  number 
of  different  appro.ximations  used  to  implement  a  numerical  solution  it  may  be 
there  is  an  advantageous  way  of  assigning  c  and  ac  so  as  to  minimize  numerical 
errors.  The  effects  of  different  values  of  c  and  a  were  not  studied  in  this  thesis. 

To  estimate  the  noise  scaling  parameter  7  we  used  the  approximation 

2 

From  the  lognormal  sample  path  7  was  estimated  to  be  .000403  and  from  the 
Rayleigh  sample  path  .000410  which  we  averaged  (.000407)  to  estimate  7. 

We  see  from  table  4.5  that  the  percentage  of  correct  decisions  are  nearly 
identical  for  the  detector  with  correct  parameter  values  and  for  the  detector 
with  estimated  parameter  values.  Under  hypothesis  1  there  is  a  slight  increase 
for  the  average  detection  time  of  11  percent  while  under  hypothesis  0  there  is 
a  decrease  in  detection  time  of  25  percent  for  the  estimated  parameter  \'alues 
compared  with  the  solution  with  correct  parameter  values. 
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4.5  Additional  observations  of  the  Zakai  de¬ 
tector 

In  addition  to  the  results  already  presented,  we  give  some  observations  made 
in  implementing  this  detector  which  one  should  be  aware  of.  We  have  already 
commented  on  the  fact  that  one  must  be  careful  in  the  choice  of  the  time  dis¬ 
cretization  step  size  At  but  one  must  also  take  care  in  the  choice  of  the  space 
(state)  discretization.  Numerical  instabilities  can  occur  for  e.xplicit  schemes  if 
the  state  discretization  step  size  is  to  small  relative  to  the  time  discetization 
size.  .Additionally,  the  computational  complexity  increases  as  the  space  dis¬ 
cretization  becomes  more  fine.  However,  if  the  discretization  is  too  coarse  then 
the  discretized  density  cannot  adequately  represent  the  actual  conditional  den¬ 
sity.  As  well  as  the  necessity  for  appropriate  choices  of  the  space  discretization 
and  time  discretization,  the  range  over  which  the  density  is  discretized  must  be 
chosen.  If  the  range  is  too  small  the  density  will  be  truncated.  If  the  range  is 
too  large  then  for  a  reasonable  number  of  points  in  the  discretization  the  den¬ 
sity  will  again  be  too  coarse.  In  figure  4.18  we  have  the  plot  of  a  signal  and 
the  conditional  expectation  of  the  signal  generated  from  a  discretization  that 
was  to  course.  One  can  see  that  the  conditional  density  could  not  adequately 
represent  the  conditional  density  and  the  conditional  expectation  of  the  signal 
tended  to  have  distinct  values  around  which  it  fluctuated.  In  figure  4.19  we  see 
the  actual  signal  and  the  conditional  expectation  of  the  signal  with  a  density 
discretized  over  too  small  a  range.  Note  how  the  conditional  expectation  of 
the  signal  cannot  track  the  signal  beyond  a  certain  point.  This  is  because  the 
discretization  range  truncated  the  density,  not  allowing  nonzero  values  for  the 
conditional  density  outside  of  the  discretized  range.  The  initialization  of 
the  underlying  Gaussian  state  densities  is  another  factor  one  should  be  aware 
of.  In  our  simulations  the  densities  are  set  to  Gaussian  distributions  using  the 
stationai-y  variance  and  mean  (whicli  is  zero).  However,  the  actual  sample  path 
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Figure  4.19:  Signal  and  conditionzd  expectation  of  the  signal  for  to  small  a 
discretization  range 
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need  not  Ijcgin  at  the  expected  value  and  it  may  take  several  iterations  before 
the  conditional  expectation  of  the  signal  yields  a  conditional  expectation  of  the 
signal  which  approximates  the  actual  signal  as  seen  in  figure  4.20.  This  can 
also  cause  initial  errors  in  the  log  likelihood  ratio  as  seen  in  figure  4.21  ,  where 
hypothesis  1  was  initialized  incorrectly.  We  see  a  general  downward  trend  for 
.01  seconds  (20  sample  points)  which  corresponds  to  the  approximate  time  it 
takes  for  the  conditional  density  to  yield  a  conditional  expectation  of  the  signal 
which  tracks  the  signal. 

4.6  The  signal  estimation  problem 

The  detection  and  estimation  problem  are  often  related.  In  this  case  we  are 
solving  the  Zakai  equation  for  the  unnormalized  conditional  density  under  two 
hypothesis  conditioned  on  observation  data.  From  the  densities  we  can  generate 
the  log  likelihood  ratio  and  implement  the  sequential  Wald  formulation  which  is 
optimal  in  the  previously  defined  sense.  However,  given  the  conditional  densities 
an  obvious  estimator  for  the  signal  is  the  conditional  expectation  of  the  signal 
conditioned  on  the  observations.  In  figure  4.22  to  4.25  we  give  the  conditional 
expectation  of  the  Rayleigh  and  lognormal  signals  under  both  hypothesis  for 
the  implicit  scheme  and  SNR=o.  We  see  that  the  Rayleigh  model,  whicli  has 
a  larger  stationary  variance,  is  better  able  to  follow  the  more  rapid  variations 
of  the  signal  while  the  lognormal  model  is  slower  and  has  more  of  a  smoothing 
effect.  We  see  that  the  solution  of  the  Zakai  equation  not  only  yields  a  sequential 
detection  scheme  but  a  signal  estimation  scheme  as  well.  We  do  not,  however, 
make  any  evaluation  of  the  perform  -3  of  this  estimator  here. 


61 


Chapter  5 
Conclusions 


We  have  jiiesented  an  "optimal”  sequential  detector  for  the  Wald  formulation 
for  diffusion  processes  for  binary  hypothesis  testing.  W^e  gave  several  numerical 
schemes  for  implementing  the  detector  for  underlying  states  of  dimension  one 
and  two  which  were  studied  to  varying  degrees.  Different  characteristics  of  the 
signal  and  discretization  scheme  were  studied  to  determine  their  effects  on  the 
detector  performance.  It  is  clear  that  the  performance  varies  significantly,  de¬ 
pending  on  many  factors  including  the  relative  decorrelation  time  of  the  signals, 
the  SNR,  the  time  discretization,  and  space  discretization  though  not  exclusively 
these  factors.  Because  of  the  significant  range  of  performance  and  the  numerous 
l)arametors  that  gave  rise  to  this  range  of  performance,  it  is  oljvious  that  one 
must  be  an  “expert”  in  order  to  make  the  detector  useful.  .Also,  the  diffusion 
model  represents  a  very  large  class  of  signals  beyond  the  lognormal  and  Rayleigh 
models  studied  here  for  which  the  detector  has  not  been  tested.  .Additionally, 
the  effect  of  many  features  of  the  signal  where  only  partially  studied.  Clearly, 
much  more  work  could  be  done  in  evaluating  the  detector’s  performance.  In 
order  to  get  a  scheme  that  is  computationally  reasonable  and  performs  well  , 
which  may  not  be  possible,  there  will  be  a  need  for  a  number  of  simulations 
with  different  parameters.  These  facts  point  to  a  need  for  automating  the  gen¬ 
eration  of  the  detector  with  “e.xpert”  supervision  for  those  not  familiar  in  the 
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theory  behind  the  Zakai  detector.  The  MACSYMA  code  is  an  initial  step  in 
automating  the  process.  The  software  system  DEsign  Laboratory  for  Process¬ 
ing  Hidden  Information  (DELPHI)  is  under  development  which  combines  an  AI 
engine,  symbolic  algebra,  and  numerical  schemes  as  an  “expert”  system  that  not 
only  generates  and  implements  the  Zakai  detector  but  other  detection  schemes 
and  model  verification  capabilities  as  well. 

We  also  noted  that  the  Zakai  detector  for  hypothesis  testing  suggests  a  signal 
estimation  scheme,  namely  the  conditional  expectation  of  the  signal  conditioned 
on  the  oliservations.  We  presented  some  results  pertaining  to  signal  estimation, 
however,  no  performance  measure  or  comparison  to  other  estimation  schemes 
was  performed. 
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Appendix  A 


This  appendix  contails  derivation  of  correlated  lognormal  and  correlated  Rayleigh 
densities. 

Consider  t  he  vector 

[.ril  .J'lo  •••  .Tim]  {•f-il  ^22  ■2^2m]  ■■  [^nl  ^n2  -^’nm] 

where  x.j  are  independent  with  respect  to  j  and  have  exponential  decorrelation 

with  respect  to  i  ie.  E  =  R  (fc)  s=  o^p^.  Then  the  covariance  matrix 

is  symmetric  with  constant  diagonals  and  m  -  1  zero  diagonals  between  each 

nonzero  one.  The  covariance  matrix  is  given  by 

-  • 

1  0  •••  0  p  0  •••  0  . 

0 

0 

p 

0 

0 
1 

Then  is  symmetric  and  tridiagonal  with  m  —  1  zero  diagonals  between  the 


Q  =  (T- 
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main  diagonal  and  the  off  diagonals  and  m  —  1  ones  on  each  end  of  the  main 
diagoiitvl.  and  is  given  b}’ 

I  1  -P  1 


1 


-p 


i  +  p- 


i+p"- 


-p 


L  1  j 

This  result  was  obtained  using  MACSYMA,  More  general  results  are  given  in 
[GL89|  Finally,  the  determinant  of  Q  is 

IQI  =  . 

So  for  the  lognormal  case  we  have 


1  p  ••• 


P 


n— 1 


Qi  =  <r- 


P 

1 


c?r'  = 


1 


a2(l  -  p) 


1  -p 
-P  I  +  P- 


l+p2  -p 

-P  1 
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and  for  the  Rayleigh  we  have 


Qr  =  (T- 


1  0  p 


,n-l 

P 

0 

1 


or^  =  -T 


1 


aHl  -  p) 


1  0  -p 

0  1  0 
—p  0  1  4- 


1  +  p-  0  -P 
0  1  0 

-p  0  1 

lQ.|=a^"(l-p2)2("-»  . 

The  distribution  for  jointly  Gaussian  random  variables  is  given  by 

Po(x)  =  - 1 - . 

(27r)t  IQI* 

Lognormal: 

Now  in  the  lognormal  case  we  have 


X  = 


(x2il  •••  (x„i 


and  the  1  to  1  transformation  yi  =  sexp(x,i).  Applying  the  Jacobian  change  of 
variables  formula  we  get 


1 


so 


nr=i  Vi 


=  1 . n)  =  Po{ln(^);i  =  l,...,u)|Jl 

s 

1 


(27r)("/2)a"(l-p2)iJifilnr»i2/i 


exp 


^  >=!*=!  L  5  J  L  S 
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where  fjjk.  is  the  jk‘^  element  of  Qi  *  which  can  be  rewritten  as 


A(z/<:?  =  i . n)  = 


exp 


-1 


2cT'^i  1  --  p^) 


Lt=l 


1=1 


n-1 


1=2 


where 


[.]  =  In(f)-/- 

=  ln(y,)  -  (ln(s)  +  p) . 


Rayleigh: 

Now  for  tlic  Rayleigh  case  we  have 

[i'll  ^12]  [i21  ^22]  •"  l^nl  in2l 
and  the  transformation  rji  =  CiJxJi  + 1?2-  Then  we  have 


PviYi  <  yi;i  =  1 


y/x;,  +r7j<yj;»=l . n 


Po{x)dx 


introducing  tlic  variable  0  and  using  the  transformation 


ill  =  l/icos(0,) 
x,2  =  yisin{9i) 


Then  the  Jacobian  is  given  by 


! 


Defining  the  vector  n  as 


[2^cos(^i)  ^sin(^i)]  [^cos(^2)  ^sin(02)]  [^cos(^„)  2^sin(^„)] 


then 


PviYi  <  yi\i  -  l.-.n)  =  L<,i.o<0<,,  Po(ll)  n  dYidOi 


70 


Assuming//  =  0  ,  expanding  the  exponent  and  using  the  identity  cos(^i )  cos(^.))+ 
sin(5i)sin(^))  =  cos(0i  —  O2)  we  get 


Pv  (V.  <  //,;  /  -  1, n)  -  (27r)("/2)a2n(l  -  /72)i"-i) 

\fl  cos(^,  -  /?,+! ) 


exp 


tldViciei 


where 


ai  =  i 


1=1  J  »=I 

-( !+/>*) 

t  =  2, n  —  1 

2<r*c»(l-/)3) 

-1 

i  =  l,n 

/2  _ 

P 

<t2c2(1  —  />2) 


Using  the  identity 


f  exp  (-4  cos(^  —  ^o)]  =  ^(-4) 

Jo<0<-2ir  2<r 


Iq  the  modified  zero  order  Bessel  function.  Then  we  get 
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Appendix  B 


Mean,  median,  power,  and  variance  of  marginal  lognormal  and  Rayleigh  densi¬ 
ties. 

Lognormal: 


7/  =  sexp(x)  ;  X  ~  N{n,a'-) 


^  1  J  N?)-M 

fiy)  =  exp  ^  - 


y/^ay 


2a2 


;  yc[0,oo 


mean 


E(y]  =  / 
Jo 


»  1  f  Hi!)-**] 
.  ^ 


2^ 


}  dy 


=  sexp 


I'*?) 


median 


r_ 

Vo  -v/i 


y/Tjray 


<'xp  <  - 


f  _  i‘<i)  1 


2a^ 


) =  /-» ■’  {-t}  =  5 


-OO  \/^a 

The  second  integral  looks  like  a  Gaussian  distribution  and  the  median  occurs  at 
the  mean  so  In  =  /t 

median  =  5exp{/i} 


power 


E(y^i  = 

Jo  y/^a  1  2a^ 


2  ^ 


y/^o 

=  s^exp{2(^-l-cr^)} 


?  dy 
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x'iinaaice 


E  [(y -.Ely])-)  =  E^  -  E|sf 

=  s'exp {2p  +  (J-)  [exp{a'}  —  l] 


Rayleigh: 


mean 


median 


power 


y  =  c\Jx\  +  x\  ;  X,  iid  N(0,a') 
/(!/)  =  ;|5«p{-5^}  ;  p^O.co] 


E|>i  =  -il}''!' 


-s 


IT 

2^^ 


r  {-2!;^} - ''‘p  {-2^} = 5 


median  =  0n(4)ccr 


=  2^^^ 


\ariaiicc 

E((y-Elyl)-)  =  Ely'^)  -  Ely]* 
(4  -  7r)c-a2 
2 
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Appendix  C 


Determination  of  stationarj'  distribution  of  the  underlying  state  equation 


dx  =  Axdt  +  DdW 


where  A  and  D  are  constants.  After  discretizing  we  have 


=  il+AAt)Xn+BAVV„  ;  AW„  iid  N(0,At) 

Denote  .4=1+  AAt  then 

x„  =  A^Xo  +  '^A’'-^-'BAWi  ;  Xo  -  iVCO.cr^) 

1=0 

Considering  the  scalar  case  which  easily  extends  to  the  vector  case  for  diagonal 
matrices  \v('  have 

E[x„]  =  0 


Denoting 


yi  =  A^-^-'BAW, 


then 


Vi  - 


and 

A^xo  ~  A(0,a2v2) 


so 


x„  -  A(0,  +  A^^^-^-'^B^At) 


75 


Now.  for  a  stationary  distribution  to  exist  must  con¬ 

verge.  So  —  1  <  ,4  <  1  or  ^  <  A  <  0.  If  we  wish  xq  to  be  initialized  to  the 
stationary  distribution  then 

.4-V2  -h  X; 


.•\..ssuming  o([nality  in  the  limit  then 


I to 


-^4^-  r  .4^’"  =  a- 

1  -  to 


With  the  prior  restriction  —1  <  -4  <  1  the  sum  converges  to  and  .4-" 


converges  to  zero,  so 


B'^At 


B^At 


l-i2  l~(H./4At)2 

So  the  diffusion  equation  has  a  stationary  distribution  of  A^(0, 
provided  ^  <  .4  <  0. 
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Appendix  D 


We  relate  the  underlying  stationary  distribution  to  underlying  state  diffusion 
parameters  for  the  stochastic  differential  equation 

(/x  =  Axdt  +  Dd\'V 


where  A  and  D  are  constants.  After  discretizing  we  have 


x„+i  =  (1  +  .4Af)x„  +  ;  AW„  iid  N(0,  At) 


Denoting  .4  =  1+  .4Af  then 


n-l 


x„  =  .4"xo  +  Y.  -4’‘"‘"'i?AlF,  ;  xo  ~  A'(0,a2)  . 


1=0 


From  .4i)p<!U(lie  C  we  have  E{x„]  =  0  so  for  the  scalar  case 


k-l 


R{,k)  =  E[x„+kX„)  =  +  Y  .VBAWn+k-i-j  Xn 

j=0 

=  •4''E(x2l 
=  . 

If  we  set  .4  =  exp|  — then  .r„  has  exponential  decorrelation  with  time  con¬ 
stant  tr  and 

p  —  exp  =  A  =  1  -  AAt 


or 


A  = 


At 
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Appendix  E 


The  following  tables  summarize  the  simulation  results  predominantly  discussed 
in  Chuprer  4.  Where  unspecified  Pp  =  P\/  =  001. 
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Hypothesis  0 

4 

iff 

%  correct 

avg.  dot. 

time  correct 

%  wrong 

avg.  dct. 

time  wrong 

avg.  det. 

time 

stat. 

loot, 

63.4 

190.7 

36.6 

1341.3 

611.4 

meams 

lot. 

31.5 

327.0 

68.5 

1129.6 

877.2 

match 

tl 

94.9 

761.1 

5.1 

994.1 

773.0 

ti/io 

100.0 

114.2 

0 

- 

114.2 

ti/100 

100.0 

29.9 

0 

- 

29.9 

loot, 

74.9 

284.8 

25.1 

1341.8 

549.9 

lot. 

36.8 

430.8 

63.2 

1301.4 

981.4 

tl 

96.8 

781.2 

3.2 

1099.6 

791.5 

ti/10 

100.0 

123.9 

0 

- 

123.9 

t,/100 

100.0 

30.2 

0 

- 

30.2 

■ 

Hypothesis  1 

B 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  dct. 

time 

stat. 

loot. 

100.0 

700.4 

0 

- 

700.4 

means 

lot. 

100.0 

780.9 

0 

- 

780.9 

match 

tl 

100.0 

936.9 

0 

- 

936.9 

ti/10 

96.4 

1317.2 

3.6 

788.6 

1298.4 

t,/100 

0 

- 

100.0 

48.8 

48.8 

loot. 

100.0 

776.4 

0 

- 

776.4 

lot. 

100.0 

915.0 

0 

- 

915.0 

tl 

100.0 

1091.2 

0 

- 

1091.2 

t,/10 

93.7 

1616.4 

6.3 

1241.2 

1592.7 

t,/100 

0 

- 

100.0 

47.8 

47.8 

Table  E.l:  Implicit  scheme,  SNR=2  Pm  =  Pp  =  .001 
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Hypothesis  0 


■ 

Hypothesis  0 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

loot/ 

63.0 

72.7 

37.0 

658.6 

289.7 

means 

lot/ 

25.5 

113.3 

74.5 

576.6 

458.5 

match 

t/ 

79.5 

487.2 

20.5 

664.7 

523.5 

t,/10 

100.0 

68.2 

0 

- 

68.2 

t//100 

100.0 

22.4 

0 

- 

22.4 

loot/ 

70.1 

75.8 

29.9 

824.3 

299.3 

lot/ 

35.6 

133.9 

64.4 

637.0 

457.9 

t/ 

83.3 

493.7 

16.7 

741.4 

535.1 

t//10 

100.0 

72.2 

0 

- 

72.0 

t//100 

100.0 

21.4 

0 

- 

21.4 

Hypothesis  1 

^  "tf 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

loot/ 

100.0 

419.7 

0 

- 

419.7 

means 

lot/ 

100.0 

468.9 

0 

- 

468.9 

matcii 

u 

100.0 

500.6 

0 

- 

500.6 

t,/10 

99.5 

393.0 

0.5 

164.4 

391.9 

ti/m 

0 

- 

100.0 

78.5 

78.5 

loot/ 

100.0 

488.5 

0 

- 

488.5 

lot/ 

100.0 

508.0 

0 

- 

508.0 

ti 

100.0 

597.0 

0 

- 

597.0 

t//10 

99.8 

456.1 

0.2 

400.0 

456.0 

t//100 

0 

- 

100.0 

73.9 

73.9 

Table  E.3:  Implicit  Scheme,  SNR=5  Pm  =  Pp  =  .001 
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■I 

Hypothesis  0 

B 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

loot/ 

63.0 

72.7 

37.0 

658.6 

289.7 

means 

lot/ 

25.8 

117.0 

74.2 

576.4 

457.9 

match 

tl 

79.3 

498.5 

20.7 

681.2 

536.4 

tt/io 

100.0 

67.0 

0 

- 

67.0 

t//100 

100.0 

22.1 

0 

- 

22.1 

loot/ 

70.1 

75.8 

29.9 

824.3 

299.3 

lot/ 

35.2 

135.5 

64.3 

637.6 

461.0 

tl 

83.9 

502.6 

16.1 

731.8 

539.4 

ti/io 

100.0 

71.2 

0 

• 

71.2 

t//100 

100.0 

21.0 

0 

- 

21.0 

Hypothesis  1 

ss  tr 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

loot/ 

100.0 

419.3 

0 

- 

419.3 

means 

lot/ 

100.0 

468.6 

0 

- 

468.6 

match 

tl 

100.0 

499.4 

0 

- 

499.4 

ti/io 

99.6 

380.0 

0.4 

153.0 

380.0 

t//100 

0 

- 

100.0 

84.7 

84.7 

loot/ 

100.0 

489.1 

0 

- 

489.1 

lot/ 

100.0 

508.1 

0 

- 

508.1 

tl 

99.9 

595.1 

0.1 

341.0 

594.9 

ti/io 

99.9 

439.5 

0.1 

309.0 

439.4 

t//100 

0 

- 

100.0 

79.2 

79.2 

Table  E.l:  Mixed  scheme,  SNR=5  Pst  =  Pf  =  .001 
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Hypothesis  0 


»  tr 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

mt, 

72.7 

16.3 

27.3 

219.2 

71.7 

means 

lot/ 

31.8 

22.7 

68.2 

209.5 

150.1 

matcli 

ti 

40.6 

94.6 

59.4 

207.0 

161.3 

tl/lO 

99.1 

28.4 

0.9 

52.0 

28.6 

ti/ 100 

100.0 

10.1 

0 

- 

10.1 

loot/ 

91.3 

15.0 

8.7 

431.0 

51.2 

lot/ 

44.6 

17.3 

55.4 

240.9 

141.1 

t/ 

50.0 

71.7 

50.0 

238.0 

154.8 

t//10 

99.5 

36.3 

0.5 

43.5 

36.4 

t//100 

100.0 

10.9 

0 

- 

10.9 

■ 

Hypothesis  1 

B 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

loot/ 

99.8 

180.6 

0.2 

584.0 

180.4 

means 

lot/ 

100.0 

199.9 

0 

- 

199.9 

match 

t/ 

100.0 

233.5 

0 

- 

233.5 

t//10 

100.0 

144.9 

0 

- 

144.9 

t//100 

92.5 

132.6 

7.5 

85.0 

129.1 

loot/ 

100.0 

206.7 

0 

- 

206.7 

lot/ 

100.0 

224.7 

0 

- 

224.7 

t/ 

100.0 

237.5 

0 

- 

237.5 

t^/io 

100.0 

181.0 

0 

- 

181.0 

t//100 

91.7 

224.8 

8.3 

248.0 

226.8 

Table  E.2:  Implicit  scheme,  SNR=20  Pm  =  Pf  =  .001 
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■ 

Hypothesis  0 

B 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

loot/ 

72.7 

16.3 

27.3 

219.2 

71.7 

means 

lot/ 

32.1 

25.6 

67.9 

217.3 

155.8 

match 

tl 

41.3 

93.9 

58.7 

204.2 

158.7 

ti/io 

99.1 

27.9 

0.9 

52.0 

28.1 

ti/lOO 

100.0 

9.9 

0 

- 

9.9 

loot/ 

99.1 

14.7 

8.1 

518.8 

55.5 

lot/ 

44.6 

17.3 

55.4 

240.9 

141.1 

tl 

50.0 

72.0 

50.0 

237.8 

154.9 

t//10 

99.5 

35.9 

0.5 

43.5 

35.9 

t//100 

100.0 

10.8 

0 

- 

10.8 

■ 

Hypothesis  1 

1 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

stat. 

loot/ 

99.8 

180.5 

0.2 

584.0 

181.3 

means 

lot/ 

100.0 

203.3 

0 

- 

203.3 

match 

tl 

100.0 

233.3 

0 

- 

233.3 

ti/io 

100.0 

136.9 

0 

- 

136.9 

tl/ 100 

95.5 

87.6 

4.5 

56.7 

86.2 

loot/ 

100.0 

202.0 

0 

- 

202.7 

lot/ 

100.0 

224.7 

0 

- 

224.7 

tl 

100.0 

274.2 

0 

- 

274.2 

ti/io 

100.0 

175.0 

0 

- 

175.0 

ti/100 

98.6 

118.7 

1.4 

37.0 

117.5 

Table  E.3:  Mixed  scheme,  SNR=20  Pm  =  Pp  —  .001 


Hypothesis  0 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

100.0 

85.5 

0 

- 

85.5 

100.0 

68.2 

0 

- 

68.2 

100.0 

49.4 

0 

- 

49.4 

100.0 

26.3 

0 

- 

26.3 

100.0 

14.9 

0 

- 

14.9 

100.0 

9.8 

0 

- 

9.8 

Hypothesis  1 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

100.0 

392.9 

0 

- 

392.9 

99.5 

393.0 

.5 

164.4 

391.9 

98.0 

371.6 

2.0 

137.5 

3G7.0 

88.8 

357.8 

11.2 

101.0 

329.0 

53.1 

281.1 

46.9 

34.7 

165.6 

41.0 

22C.1 

59.0 

22.1 

105.7 

Table  E.7;  Implicit  scheme,  SNR=5,  varying  the  probability  of  a  miss 
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Pf 

Pm 

Hypothesis  0 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

.0001 

.001 

100.0 

69.9 

0 

- 

69.9 

.001 

.001 

100.0 

68.2 

0 

- 

68.2 

.01 

.001 

100.0 

66.1 

0 

- 

66.1 

.1 

.001 

08.1 

66.0 

1.9 

51.4 

65.7 

.5 

.001 

G3.9 

47.8 

36.1 

10.8 

34.5 

.75 

.001 

36.4 

38.8 

63.6 

3.6 

16.4 

Pf 

Pm 

Hypothesis  1 

%  correct 

avg.  det. 

time  correct 

%  wrong 

avg.  det. 

time  wrong 

avg.  det. 

time 

.0001 

.001 

99.8 

506.2 

.2 

332.0 

505.8 

.001 

.001 

99.5 

393.0 

.5 

164.4 

391.9 

.01 

.001 

99.7 

264.4 

.3 

109.3 

263.9 

.1 

.001 

99.9 

129.4 

.1 

222  0 

129.5 

.5 

.001 

99.9 

26.2 

.1 

370.0 

26.5 

.75 

.001 

99.8 

6.8 

.2 

86.5 

7.0 

Table  E.S:  Implicit  scheme,  SNR=5,  vaning  the  probability  of  false  alarm 
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